%This is the Matlab code used to create Table 1 in "Why experimenters might not always want to randomize, and what they could do instead"

clear;
n=4;
ind = (0:2^n-1)';
D = de2bi(ind)
X=(0:n-1)';
nd=sum(D,2);
variance=1./nd + 1./(n-nd)

f1=X';
f2=-X'.^2;
bias1=D*f1'./nd - (1-D)*f1'./(n-nd);
bias2=D*f2'./nd - (1-D)*f2'./(n-nd);
MSE1=bias1.^2+variance;
MSE2=bias2.^2+variance;

[MSEsort, order]=sortrows([variance MSE1])

designs=[ones(1,16);
         [ones(1,14) zeros(1,2)];
         [ones(1,6) zeros(1,10)];
         [ones(1,4) zeros(1,12)];
         [ones(1,2) zeros(1,14)]];
designsnorm=designs./(sum(designs,2)*ones(1,16))


table=[D(order,:) designs' bias1(order) variance(order) MSE1(order) bias2(order) variance(order) MSE2(order) zeros(2^n,1)]
     
MSE1designs=designsnorm(:,1:n^2-2)*MSE1(order(1:n^2-2))
MSE2designs=designsnorm(:,1:n^2-2)*MSE2(order(1:n^2-2))
MSE1designs(1)=Inf;
MSE2designs(1)=Inf;

%%
%calculate our criterion function manually
Xnorm=(X-mean(X))/std(X);

C=zeros(2*n);
for i=1:n
    for j=1:n
        C(i,j)=10*exp(-(X(i,:) - X(j,:))* (X(i,:) - X(j,:))'/10);
    end
end
C(n+1:2*n,n+1:2*n)= C(1:n, 1:n);
C(1:n,n+1:2*n)= C(1:n, 1:n)*exp(-1/10);
C(n+1:2*n,1:n)= C(1:n, 1:n)*exp(-1/10);


 

for i=1:2^n
    d=D(order(i),:);
    n1=sum(d);
    n0=n-n1;
    w=[(1/n1)*d - (1/n)*ones(1,n)  -(1/n0)*(ones(1,n)-d) + (1/n)*ones(1,n)];
    
    ebias2=w*C*w';
    table(i, 16)=variance(order(i)) + ebias2;
end
variance(order(:))
table(:,16)


%%
FID = fopen('designcomparisonSimple.tex', 'w');
fprintf(FID, '\\begin{tabular}{cccc|ccccc|ccc|ccc|c}\\hline \\hline \n');
fprintf(FID, '\\multicolumn{4}{c}{assignment} & \\multicolumn{5}{c}{designs}  & \\multicolumn{3}{c}{model 1}  & \\multicolumn{3}{c}{model 2}  & EMSE   \\\\ \n');
fprintf(FID, ' $d_1$ & $d_2$ & $d_3$ & $d_4$  &  1 & 2 & 3 & 4 & 5 & bias & variance & MSE  & bias & variance & MSE &     \\\\ \n');

fprintf(FID, '\\hline \n');

for d=1:n^2
    fprintf(FID, '%i & %i & %i & %i  & %i & %i & %i & %i & %i& %8.1f &%8.1f &%8.1f  &%8.1f &%8.1f &%8.1f  &%8.2f\\\\ \n ', ...
    table(d, 1),  table(d, 2),  table(d, 3),  table(d, 4), table(d, 5),  table(d, 6),  table(d, 7),  table(d, 8), ...
    table(d, 9),  table(d, 10),  table(d, 11),  table(d, 12), table(d, 13),  table(d, 14),  table(d, 15),  table(d, 16));

    if (d==6)
        fprintf(FID, '\\hline \n');
    end
end
fprintf(FID, '\\hline \n');
fprintf(FID, '\\multicolumn{4}{c|}{MSE model 1:} & %8.1f &%8.1f &%8.1f &%8.1f &%8.1f & & & & & & \\\\ \n ', ...
        MSE1designs(1), MSE1designs(2), MSE1designs(3), MSE1designs(4), MSE1designs(5));
fprintf(FID, '\\multicolumn{4}{c|}{MSE model 2:} & %8.1f &%8.1f &%8.1f &%8.1f &%8.1f & & & & & & \\\\ \n ', ...
        MSE2designs(1), MSE2designs(2), MSE2designs(3), MSE2designs(4), MSE2designs(5));    
    

fprintf(FID, '\\hline \\hline \n');
fprintf(FID, '\\end{tabular}\n');
fclose(FID);
